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ABSTRACT 

Electron-positron pairs may be produced near accreting black holes by a va- 
riety of physical processes, and the resulting pair plasma may be accelerated and 
coUimated into a relativistic jet. Here we use a self-consistent dynamical and ra- 
diative model to investigate pair production by 77 collisions in weakly radiative 
accretion flows around a black hole of mass M and accretion rate M. Our flow 
model is drawn from general relativistic magnetohydrodynamic simulations, and 
our radiation field is computed by a Monte Carlo transport scheme assuming 
the electron distribution function is thermal. We argue that the pair production 
rate scales as r~^ . We confirm this numerically and calibrate the scaling 
relation. This relation is self-consistent in a wedge in M, M parameter space. If 
M is too low the implied pair density over the poles of the black hole is below the 
Goldreich- Julian density and 77 pair production is relatively unimportant; if M 
is too high the models are radiatively efficient. We also argue that for a power- 
law spectrum the pair production rate should scale with the observables Lx = 
X-ray luminosity and M as L\M~^. We confirm this numerically and argue that 
this relation likely holds even for radiatively efficient flows. The pair production 
rates are sensitive to black hole spin and to the ion-electron temperature ratio 
which are flxed in this exploratory calculation. We flnish with a brief discussion 
of the implications for Sgr A* and M87. 

Subject headings: accretion, accretion disks — black hole physics — MHD — 
radiative transfer — Galaxy: center 



1. Introduction 



Models of zero-obliquity black hole accretion — in which the accretion flow angular mo- 
mentum is parallel to the black hole spin — typically exhibit a low density "funnel" over 
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the poles of the black hole. The funnel is empty because the funnel plasma is free to fall 
into the hole or be ejected to large radius. Magnetic fields do not prevent this: in mag- 
netohydrodynamic (MHD) simulations of radiatively inefficient accretion flows (RIAFs) the 
funn el magnetic field typically runs in a smooth spiral from the event horizon to large ra- 



dius ( iDe Villiers et a 



2006 



Beckwith et al. 



20031 : iMcKinnev &: Gammiell2004l : lKomissarovll2005l : iHawlev fc Krolik 
20081 ). Because the field lines do not leave the funnel there is no way 



for the disk plasma to resupply the funnel plasma. 

What process, then, populates the funnel with plasma? And what controls the tem- 
perature (or distribution function) of the funnel plasma? These questions bear directly on 
two interesting problems in black hole jet theory: are jets made of pairs or an electron-ion 
plasma? And which is more luminous: the base of the jet or the accretion fiow? The purpose 
of this paper is to investigate these questions in the specific context of hot, underluminous 
accretion fiows where nearly ah initio models are computationally feasible. 

There are several pair creation processes that might populate the funnel with plasma. 
Plasma close to the event horizon in a RIAF is relativistically hot, and thus can form 
electron-positron pairs through particle-particle (ee, ep), particle- photon (e7, ^7), or 
photon-photon collisions (77). The cross section near the energy threshold is largest 
for 77 interactions, which have a cross section ~ o"t =, the Thomson cross section. In 
the funnel the photon density vastly exceeds the particle density, so 77 collisions dom inate 



production (IStepney fc GuilbertI Il983l . iPhinneyl llQSSl. iPhinnev 



1995 



prod uction by these processes is discussed in e.g. Kusunose &: Mineshigel 



Krolik 



19991 ). Pair 



19961 and 



Esin 



1999 in the context of Advection Dominated Accretion Flows (ADAFs, Narayan fc Yi 19941 ). 



These works, however, focus on the energetic role of pairs in ADAF disks rather than the 
population and dynamics of pairs in the funnel. 



Other processes are important when the density is below the iGoldreich fc Julian! (119691 ) 
charge density. Then the plasma can have E ■ B 7^ 0, and the electric field can directly 
accelerate particles to high Lorentz factors. The energetic particles Compton upscatter 
background photons that colli de with other background photons and produce a shower of 
pairs in a pair-photon casc ade f Blandford fc Znaiek 1977: |Phinnevlll983l : iBeskin et al.lll992 
Hirotani fc Okamotd Il998l . and recently IVincent fc Lebohed I2OIOI ) . 



In this paper we model production of an plasma by photon-photon collisions in the 
funnel above a hot, underluminous accretion disk. At low accretion rates M (< Merit ~ 
10~^L£;rfrf/(0.1c^), where Lem = Eddington luminosity) the disk cools on a timescale longer 
than the accretion timescale; it is a RIAF. In this regime the radiative and dynamical evolu- 
tion are decoupled and it is practical to treat both on a nearly ab initio basis. Throughout 
the range of M we consider the funnel pair plasma is tenuous enough that annhilation is 
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negligible, so pair production will be balanced by advective losses such as accretion into the 
black hole or loss in a wind. 



We draw our RIAF model from two and three dimensional gen eral relativistic 



netohydrodynami cs simulations (GRMHD, using the HARM code, iGammie et al.l 12003 



mag 



Noble et al.ll2009l ) of an accreting, magnetized torus with zero cooling. The ra diation field is 



calcu lated as a post-processing step using a Monte Carlo method (grmonty, iDolence et al. 



20091 ). and pair production rates are estimated from snapshots of the radiation field using a 



procedure described in detail below. 

This paper is organized as follows. In § [2] we describe the basic model for accretion 
flow dynamics and radiative transfer. In § [3] we write down the pair production model and 
present a test problem for our Monte Carlo scheme. Scaling formulas are presented in § |H In 
§ E] we show results for a range of black hole masses and accretion rates. We briefly discuss 
implications for Sgr A* and M87 in § [6] and summarize in § [71 



2. Accretion flow model 

We use a numerical model for the accretion flow and for the radiation field; together 
these nearly ab initio models form a numerical laboratory for investigating physical processes 
near a black hole in a self-consistent way. 



2.1. Dynamical model 

We use a relativistic MHD model for the accret ing plasma (see e.g., 



Gammie et al. 



20031 ) . The initial condition is an equilibrium torus (jFishbone fc Moncriej Il976l ) in orbit 
around a Kerr black hole with a^, = 0.94, where a^GM'^/c is the hole angular momentum. 
The torus is seeded with poloidal, concentric loops of weak magnetic field that are parallel 
to density contours. Small perturbations are added to the internal energy and this seeds 
the magnetorotational instability, which leads to the development of MHD turbulence in the 
disk and accretion onto the central black hole. The model extends from slightly inside the 
event horizon to r = AOGM/c^. 

We solve the evolution equations until a quasi-equilibrium accretion flow is established, 
meaning that the mean structure of the flow is not evolving on the dynamical timescale. 
Our (untested) hypothesis is that at r < IbGM/c^ the model accurately represent the inner 
portions of a relaxed accretion flow extending over many decades in radius. 
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A few of the physical assumptions in the GRMHD model are worth stating explicitly. 
The equation of state is 

P = {lad - 1)m (1) 

where 7ad = 13/9 (appropriate for ion temperature Tj < nipC^/k = 1.1 x 10^^ K and electron 
temperature Tg > meC^/k = 5.9 x 10^ K), p = pressure, and u = internal energy density. 
Particle number is also conserved: 



0, 



(2) 



where po = rest-mass density and = four-velocity, in the dynamical evolution. That is, 
pair production is not included in the dynamical model. The model is therefore consistent 
only if pair creation is weak enough not to alter the flow dynamics or energetics. 



We evolve the GRMHD equations using the harm code ( IGammie et al.ll2003l : iNoble et al. 



20091 ). harm is a conservative scheme that evolves the total energy rather than internal energy 
of the flow. The M HD equation integrati on is performed on a uniform grid in modified Kerr- 
Schild coordinates (IGammie et al.ll2003l ). The coordinates are logarithmic in Kerr-Schild 
radius r and nonuniform in Kerr-Schild colatitude 6 (Boyer-Lindquist and Kerr-Schild r and 
6 are identical), with zones concentrated toward the midplane of the accretion disk. The inner 
and outer radial boundaries use outflow boundary conditions. The axisymmetric models use 
a 256 X 256 grid, and the single 3D run uses a 192 x 192 x 128 gr id. For detai ls of the 



num erical method, the init i al setu p and the flow evolution in 2D see IGammie et al.l (j2003[ ) 
and iMcKinney fc Gammid (j2004f ). A snapshot of the density, temperature, and magnetic 
field strength from one of our runs is shown in Figure [TJ 



2.2. Radiative model 



Our radiative model is identical to that applied by lMoscibrodzka et al.l ( 120091 ) to Sgr A*, 
although here we restrict attention to a thermal plasma with Te = (except in one case 
noted below). Synchrotron emission and absorption are included, as is Compton scattering. 

Bremsstrahlung is not important in the inner parts of the accretion flow. For a thermal 
plasma with Gg = kT^/ {meC^) > 1 the ratio (synchrotron / bremsstrahlung) cooling ~ 
Gg/(a/3), where a = fine structure constant and /3 = Sirp/B"^. At the radii of interest 
here 6e ~ 1 — 10^, and /3 ~ 10, so in an energetic sense synchrotron dominates the direct 
production of photons. 

Synchrotron emission occurs at a characteristic frequency Us ~ {eB / {2nmec))Ql which 



- 5 - 



is <C rrieC^/h for any astrophysically reasonable combination of M and M 0. Potentially 
pair-producing photons must therefore be produced by Compton scattering. 

We compute the r adiation field using the general relativistic Monte Carlo radiative 



transfer code grmonty (IDolence et al.ll2009l ). The radiation field is represented by photon 
packets (photon rays or "superphotons" ) . Each superphoton is characterized by a weight 
w = number of physical photons/superphoton, and a wave four- vector k'^. Superphotons 
are produced by sampling the emissivity. The wavevector is transported according to the 
geodesic equation. Along a geodesic w is decremented to account for synchrotron absorption. 
Compton scattering is incorporated by sampling scattering events. When a superphoton 
scatters it is divided into a scattered piece with new wavevector k''^ and new weight w', and 
an unscattered piece along the original wavevector with weight w — w'. The distribution of 
scattered k''^ is consistent with the full Klein-Nishina differential cross section. 

We use a "fast light" approximation in treating the radiative transfer. The data from 
a single time slice t„ (e.g. po(tn,x^,x'^,x^)) is used to calculate the emergent radiation field 
as if the data, and therefore photon field, were time-independen t. We have checked the 



fast light model against a time-dependent radiative transfer model (IDolence et al.ll2010[ ) and 
verified that this approximation does not introduce significant errors. 



2.3. Model scaling 

The properties of the accretion flow model are independent of the absolute value of 
the density (provided the magnetic field strength is scaled appropriately), but the radiative 
model is not. To scale the model we specify the length unit 

C . (4) 



For the synchrotron emissivity we use the approximate expression of iLeung et Zl |2010l ) 



= ^^|^(xV2 + 2"/i2xV6f exp(-XV3) 

3cK2{Qe ) 

where X = vjvs^ = 2/9(ei?/27rmec)0g sin0 is the synchrotron frequency, d is an angle between the 
magnetic field vector and emitted photon, and is a modified Besscl function of the second kind. The 
fractional error for this approximate formula is smaller than 1% for Og > 1 (where most of the emission 
occurs) and increases to 10% and more at low frequencies for 8e < 1 (where there is very little emission). 
The synchrotron emissivity function peaks a% v k, %Vs. 
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time unit 



r 



GM 



(5) 



and mass unit Ai, which is proportional to the mass accretion rate. M does not set a mass 
scale because it appears only in the combination GM. Since Ai <^ M the accretion flow 
does not affect the gravitational field. 

Given M and Ai the radiative transfer calculation is well posed. Typically M can be 
estimated directly from observations, while Ai is varied until the model submillimeter flux 
matches the observed flux. 



2.4. Model limitations 

An important limitation of our model is that it treats accreting plasma as a nonradiating 
ideal fluid. This implies that electrons and ions have an isotropic, thermal distribution 
function. The potentially important effects of pressure anisotropy and conduction (e.g. 
Sharma et al.ll2006l . I Johnson &: QuataertI 120071 ) are therefore neglected, as are the radiative 
effects of a nonthermal component in the electron distribution function. 



Cooling is also neglected. This is a good approximation in low accretion rate systems 
like Sgr A*, but a poor approximation in higher accretion rate systems like M87. If one 
were to turn on cooling but hold the synchrotron flux flxed the density and magnetic fleld 
strength (i.e. the mass unit Ai) would increase. 



3. Pair production 



3.1. Basic equations 

For a population of photons with distribution function dN^/d^xd^k (here d^k = dkidk2dks 
and 1, 2, 3 are the spatial coordinates) the invariant pair production rate per unit volume is 



dN. 



± 



'—g d^xdt 2 



d^k d^k' dN^ dN, 



7 J2 



1 fX'Y'Y C 



^k* ^k't d^xd^k d^xd^k' l^^^l 



(6) 



where g is the determinant of g^^u, and the factor of 1/2 prevents double-counting. Here 
is the cross section for 7 + 7 — )■ e"*" + e~: 



a. 



77 



;[CM] 



(2e[cM]^ + 2e[cM]^ — 1) cosh ^ e[cM] — e[CM](e[CM]^ + 1)\/ e[CM]^ — 1 (7) 
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^[CM] = -UcMf,k^ = -U[CM],j,k'^ 



(Breit & Wheeler 1936), 

is the energy of either photon in the center of momentum ([CM]) frame of the two photons, 
and u\cM] is the four-velocity of the [CM] frame. 

Equation E] is coordinate invariant since ^—gdP'xdt is invariant, the distribution function 
is invariant (because d^xd^k is invariant), epM] is a scalar, the cross section is invariant, and 
d^k/ y/—gk^ is invariant. It reduces to the correct rate (cf. eq. 12.7 of Landau & Lifshitz, 
Classical Theory of Fields) in Minkowski space, and is therefore the correct general expression 
for the pair production rate. Because n± itself is invariant it also describes the pair creation 
rate in the fluid frame. 

We will need the rate of four-momentum transfer from the radiation field to the plasma 
via pair creation: 

1 dP^ _ ^1 f d^k d^k' dN^ dN^ 2 
-^d^xdt~ 2j ^k^^k'^dHd^kd^xd^k'^ + )e[cA/]^77C- (9) 

Here A is a constant that makes the equation dimensionally correct. 



3.2. Monte Carlo estimate of pair creation rate 

We estimate the integrals (j6]) and ([9]) using a Monte Carlo scheme. Given a sample of 
photons on a time slice t within a small three- volume A^x, a naive estimate is 

1 dN± 1 sr-^ f \ f "^j \ ^ 1 2 ^n^ 

^d^xdt 2 ^ I A% J I A% J T^T^'I^A^]^^^'^- ^^"^ 

where i and j label superphotons. 

If there are A^^ superphotons in A^x then there are 0{Ng) pairs of superphotons and the 
computational cost of ( ITU]) is 0{N'^). One might hope that the error would scale as 1/ ^/N^ 
because there are 0{Ng) pairs, but this is wrong. There are only A^^ independent samples 
and so the error scales as 1/ y/W^. 

We obtain an estimate with accuracy that is the same order as flTUl) at 0{Ns) cost by 
selecting an unbiased sample of A'^s pairs of superphotons and evaluating 

^d^xdt ^ 22 ^ \AH) \AH) J^k\ ^kf^""''^ ^ ^ 

i,j V ^ J V ^ 



An identical procedure is used to evaluate G^. 
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3.3. Test problem 

Does our Monte Carlo procedure accurately estimate the pair production rate? As a 
check we evaluate pair production rates near two point sources of £" = 4meC^ photons. The 
calculation is done in Minkowski space and Cartesian coordinates (so a/— g = 1), and the 
optical depth to pair creation is assumed small. At each point we compare the analytic and 
numerical result. 

The expected pair production rate is given by Equation [61 The energy of two colliding 
photons in their center-of-momentum frame is a function of the cosine /i of the angle between 
the rays from the two sources: e^-^^j = (1/2) (1 — fi)k'^k"^. The photon momentum space 
distribution is a 5 function. The number density of photons dN^/d^x = N^/^Anr'^c) at 
distance r from the source, where each source produces photons at rate A'^. 

Figure [2] shows a 2D map of the numerically evaluated pair production rate in the plane 
of the two sources. Figure [3] shows the analytic and numerical pair production rates along 
the black contour shown in Figure [2] (upper panel), and their difference (lower panel). The 
error in the numerical rate is oc A^s , as demonstrated in Figure HI where A^, is the number 
of photon packets emitted by each source. Evidently the Monte Carlo method produces an 
unbiased, convergent estimate of the pair production rate. 



4. RIAF scaling laws 

In this section we derive scaling relations for the pair production rate in two cases: (1) 
M and M are known and the flow is radiatively inefficient; (2) the spectrum uL^, and M are 
known. In case (1) we can numerically evaluate the pair production rate self-consistently and 
check it against the scaling relation. In case (2) we can do the same, but we also obtain a 
method for estimating pair production rates from observations that may also apply to flows 
that are radiatively efficient. 

First consider the pair production rate density at a single point in the flow where the 
plasma-frame photon spectrum is a power-law with high energy cutoff at e = E/{meC^) = 

^max 3> 1. 

Hn. n n I 



dE rrieC^ 

We evaluated the pair production rate density numerically for this energy distribution. A 
flt to the result over — 3 < a < 2 and 10 < e^ax < 160 gives 



- ^^"/^(^+C4irin(^) (13) 



riQaTC 16 3 
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(Zdziarski [1985] gives a similar expression in the emax ^ 1 limit). At worst the fit is ~ 2 
too small for a ^ and e,nax = 160. For a < —2, which is typical of our models, the relative 
error is smaller than 60%. 

For a > {dlnuLiy/dlnu > 2) pair production is dominated by photons with e ~ emax, 
and h± is therefore sensitive to emax- For a < 0, pair production is dominated by pairs 
with e ~ 1 in the center-of-momentum frame. In this case there is an equal contribution 
from each logarithmic interval in energy in the plasma frame, and the pair production rate 
density is therefore weakly (logarithmically) dependent on e^nax- 

Our models have a < —2 so h± is insensitive to emax- Therefore the effective number 
density of pair producing photons is uq ~ L5i2/(47r£^meC^), where L512 = z/L,^(512keV) and 
n± ~ ?7,q(Ttc. Then 

where / is a dimensionless function of Kerr-Schild radius r and colatitude 6 = cos~^ fi. 

What do we expect for the spatial distribution of pair production /? The pair-producing 
photons are made by upscattering synchrotron photons in a ring of hot gas near the innermost 
stable circular orbit (ISCO) . Away from this ring the density of photons will fall off as ~ 1/r^. 
The pair production rate also depends on the angle ip between photon trajectories in the 
coordinate frame through the geometrical factor e[cM]^/^i^2 oc 1 ~ costp. At large r < 1 
so 1 — cosip oc 1/r^. Then n± ~ r^^. Compton upscattered photons are also beamed into 
the plane of the disk by the relativistic orbital motion. If the intensity of photons to be 
scattered is nearly independent of 6 then the pair production rate density should fall off 
away from the midplane as the density of upscattering electrons ~ exp(— /i^/(2cr^) where 
(Tp ~ 0.3. Gathering these estimates together we expect / ~ exp(— /i^/(2crpm^))/r^ where 



4.1. Scalings with model parameters 

Now suppose we know the mass M = mgMg (Mg = IO^Mq) and the accretion rate M = 
rhMEdd, where MEdd = LEdd/ (erefC^) and eref = 0.1 is a reference accretion efficiency. We 
assume that photons are produced in a low frequency synchrotron peak and then scattered 
to ~ 512 keV by Compton scatterings, where Ugc is 1 or 2. 

For a plasma that is optically thin to synchrotron absorption at peak, the total num- 
ber of synchrotron photons at the peak frequency produced per unit time is N^, ^ 
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^T^Vpeakjvp^ak^^ / ^^^pe^ak), wheie ju^^^^. IS the syiichrotroii emissivityB The number density of 
synchrotron photons is then n^^^^^ ^ Ni^^^^^/{47rC'^c). 

A fraction t"'^" of the peak photons are upscattered to 512 keV, where r = axUeC is the 
Thomson depth of the plasma, so ^512 = f^up^ak'''"""' ■ The mean number of Compton scatter- 
ings is Use = ^og{meC^ / hupeak) / ^ogA, whcre A ^ 166^ is the photon energy enhancement in 
single scattering by a relativistic electron, so 

risc^ ai + a2\og—. (15) 
m 

We determine Oi and 02 numerically but for a reference model with m = 10~® and mg = 
4.5 X 10~^ (Sgr A*), the average value of n^c ~ 1 — 2. 

Assuming M ~ AirpcC, the magnetic pressure is comparable to the gas pressure and 
both are ~ pc^, the plasma density, magnetic field strength, and plasma temperature (close 
to the virial temperature) scale as 

c \ m 



ne^ 7771— — (16) 



B 



2 



Svr e, 

and 



1 / rripC^ \ ( ^\ 



(17) 



Be ^ -. 

30 nip 



The mean emission Gg corresponds to the mean value near the ISCO, and therefore increases 
with a^. Combining, 

^± - ^ f 4f ) 4"^^^'^"/ - ^'^"^^ /(7'/^)' (19) 

where ^ is a dimensionless constant to be determined numerically, tq = classical electron 
radius, and «/ = the fine structure constant. From now on unless stated otherwise we will 
set ere/ = 0.1 and the mean number of scatterings n^c =1.5 (numerical results, below, show 
1.4 < rise < 1-6 for relevant M, M). Then 

n± ~ 9 X 10^^ A mg ^ /(^, /i). (20) 
To estimate the jet kinetic luminosity we need the pair production rate: 

J3„ 



iV± = / y=gf^ a;n± (21) 



"^iv^^^^ ~ 8\/2e^neB/(27TOeC^), see lLeung et al.ll2010 
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where Vhor is the horizon radius. Then 

N± ~ h±C^ ~ lO^Mm^+^^^^m^ s"^ (22) 

Only pairs made inside the funnel and at r > r^^ = stagnation radius can escape to large 
radius (they are "free pairs"); those made at smaller radius or inside the accretion flow are 
advected into the hole. The free pair fraction is therefore a small multiple of f l22|) . 

Evidently the pair production rate density is sensitive to the mass accretion rate, n± ~ 
m^. This steep dependence shuts off pair production at low accretion rates, making it difficult 
for low m systems like Sgr A* to populate their funnel with pairs. Below we will show that 
the implied funnel pair density for Sgr A* falls below the Goldreich- Julian charge density 
(see g53]). 



4.2. Scalings with observables 

Assume that from observations we know M, the X-ray luminosity Lx = IxLq (assuming 
isotropic emission), and the spectral index a = dlog{i'L^) / dloguu Self-consistent models 
then permit us to calibrate the relation between these quantities and the pair production 
rate density. Since this relation depends only on the distribution of pair-producing photons 
within the source, it seems likely that it can be applied to sources with M > Merit, in which 
cooling is important. 

If the spectrum is power-law from the keV to MeV energy, 

L,,2{Lx) ^ Lxe'-''", (23) 



where S is a constant to be determined numerically and {c'axL'^Q/TnlG'^M^) 10^*^ cm~'^s~^. 
This assumes that the observed spectrum and the plasma-frame spectrum near the black hole 
are identical. We checked the plasma-frame spectrum and found it to be a slightly blueshifted 
version of the observed spectrum; the blueshifting does not change the scaling relation. 

The pair production rate is 

^ h^C' ^ B ( f;^; 3 ) ll e^-^^" m^-i (25) 



is used to extrapolate the spectrum from keV to MeV energies. It can be evaluated from X-ray data 
but this can be inaccurate due to localized spectral features. It can be evaluated more accurately from 
millimeter /X-ray colors. 
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where {arL'^/TnlGM^c') = 10^^ s ^. The dependence on black hole mass changes between 
Equations f l2^ and (125|) because C oc mg. 



5. Pair production in RIAF - numerical results 

We now evaluate the pair production rate numerically, check whether it matches the 
expected scaling laws, and evaluate /(r, /i). To do this, we have run simulations with a range 
of M and M, assuming that the models have equal ion and electron temperatures, = Ti. 
A list of model parameters is given in Table [TJ 



5.1. Pair creation rate 

5.1.1. Dependence on model parameters: rh, m 

The n± in models A through H (see Table [T]) is well fit by 

n±(r, /x) = 3 X lO'^m^^+'^-mg ' x (L)-6e-MV(2-i) cm-^ s'S (26) 

or ^ ~ 3 in equation ( l20l) . The constant in Equation [26] is derived from models with 
Tj/Te = 1. The constant is sensitive to Ti/Tc, for T-jTc = 3 it is 10~^ times smaller. As 
expected, h± ~ at large r; surprisingly, however, this is also good fit at all r. 

The pair production scale height a± ~ 0.3 independent of rh, m. This is nearly identical 
to (Tp, the plasma scale height. Notice that a± also controls fjet the fraction produced inside 
the funnel. The funnel wall is at 

^ > " r + AGM/c^ • 

and fjet ~ 10%. Figure [5] shows a 2D contour map of h± corresponding to model C and 
Equation (l26l) . and a contour marking the approximate funnel boundary. The grid averaged 
fractional difference between time averaged MHD models A through H and Equation (l26l) 
is < 60%. Since h± is a steeply declining function of /x^, almost all free pairs are made near 
the funnel walls. 

The pair production rate is well fit by 

iV± = 4 X 10®° m^+2n.c ^2 g-i ^28) 

where a fit gives 

nsc = l + 0.03 ln(m8/m). (29) 
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Figure [6] compares the time-averaged numerical N± to Equation [281 

5.1.2. Dependence on Ix, ct, m 

The self-consistent radiative model enables us to calculate the emergent spectrum, from 
which we can measure a 2 — lOkeV luminosity Ix and a spectral slope a. The pair production 
rate density can be measured in the same models. The numerical results are well fit by 

or i3 ^ 1 in Equation ( l24l) . The fractional error of the fit is < 50% for time averaged models 
A-L. 

The pair creation rate is well fit by 

iV± = 5xlO='0/^e9-26"m8-i s"^ (31) 

Figure [7] compares N± to the semianalytic formula given by Equation [3l] for different snap- 
shots of the simulations with different mass accretion rates (models A-C) and black hole 
masses (models F-H). The semianalytic and numerical results agree well, and the scaling 
constants are close to those estimated in ^ 

Although Equations ( 1301) and ( 1311) are strictly valid only for radiatively inefficient flows 
with rh < rhcrit, they depend mainly on the geometry of the radiation field and not on 
the radiative efficiency of the flow. We speculate that they provide a good estimate of the 
pair production rate even in more efficient systems, if a± is set to the scale height of the 
Comptonizing corona and the spectrum extends to sufficiently high energy. 

5.2. Pair power and electromagnetic luminosity of a funnel 

We define the funnel pair creation "luminosity" 

L± = /,,jiV±2meC%t (32) 
where Tjet is the jet bulk Lorentz factor at large r (assuming cold flow). Then 

L± ^ 6 X lO^Vjrf m3+2"==m2r,^, ergs s'^ (33) 

and 

L± ~ 10^V,et/xe'-''"m8-ir,etergss-^ (34) 
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It is interesting to compare this with the Blandford-Znajek (BZ), or electromagnetic, lumi- 
nosity of the funnel 

Lbz = 27r [ dd^Tl (35) 

where Tl = iP'u^ut — h^ht is the electromagnetic part of the stress-energy tensor and is 
computed directly from the simulation data. The BZ luminosity is well fit by 



8 X 10^^(1 - ^l-alfmm^ ergss"^ (36) 



The scaling with a* is taken from Equation (61) of iMcKinney fc Gammid l2004l . which is a 
fit to numerical data. 

For mass comparable to that of Sgr A* and a* = 0.94, Lbz > Lj-fTjet for m < rhcrit ~ 
10~^ (for rhcrit see § 16.11) . At low accretion rates the BZ luminosity completely dominates 
the pair luminosity, because the pair luminosity is such a steep function of accretion rate. 
Because Lj- ~ mg while Lbz ^ ''^S; the rh at which Lbz ~ is higher for lower mass black 
holes. L±/Lbz ratio is shown in the Tabled] For reasonable Tjet the funnel luminosity is 
therefore electromagnetically dominated for radiatively inefficient flows. 



5.3. Energy-momentum deposition 

Some of the pairs created in the funnel will escape to large radius and some will fall into 
the black hole. In an MHD model, the escaping fraction and asymptotic Lorentz factor will 
depend on the run of pair creation rate with radius, the magnetic field structure, the energy 
density of the pair plasma, and the pair-creation four- force C^. 

Based on numerical calculations the spatial distribution of is well fit, with x = r/C, 

by 

GU(n.) = G»^«^^iMSr (37) 



M X M/{C^V 

G,,,,(x, /ij - G ^ - j^i^c?r^^ ^^^^ 

where Ai,T,C, and rij- are given in cgs units. The four-force vector components are given 
in a Kerr-Schild coordinate basis and in code units; we divide in cgs units by the unit of 
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the four- force density Ai/C'^T'^. All components of the four- force depend more steeply on 
radius than h±. The radial component of the four- force is positive at large radius, zero at 
Xst ~ 3;/5(7o(l + /^^/2), and negative at small radius. The sign of changes at the equatorial 
plane, as it should. 

Pairs are created in the funnel with an initial distribution function, which is immedi- 
ately isotropized with respect to rotation around the magnetic field. Later evolution of the 
distribution function depends on ill-understood relaxation processes; the pairs may not relax. 
Whether or not relaxation occurs the initial mean energy of the particles is of interest. So: 
what is (iri±/(i log 7e±[i?_F] (= the energy distribution of newly created pairs), where 7e±[FF] 
is the Lorentz factor of new particles measured in the fluid frame [FF] ? 

Four-momentum is conserved in pair creation, so the average Lorentz factor of the new 
leptons in the fluid frame is 

leHFF] = -\u,{k^ + k'^) (41) 

where is four-velocity of the background plasma and we assume that is in units of 
rrieC? . Figure [H] shows dh±/ d\og'~^(,±[FF] in model C at r = rh.Tisco.'^GM / (averaged over 
20 radial zones from 6* = — 10 deg; here = event horizon radius). The distribution is fiat 
and cuts off at 7e±[FF]maz ~ 100 at r = r^^co- 

If the thermalization timescale is short (which, given the low density and likely high 
temperature of the plasma, seems unlikely to us), then the rate of internal energy injection 
due to creation is m = e-t, where the kinetic energy density injection rate, in the plasma 
frame, is0 

/diTi 
dle^lFF] (7e±[FF] " l) "^eC^ 3 — (42) 

The corresponding temperature of the newly injected pairs is 6e± = (l/3)e±/(n±meC^) ~ 
(20, 10,3) at r = (r/^, rjgco, 5 GM/c^). This temperature is comparable to the thermal back- 
ground plasma and is sub-virial. Pairs are not born hot. 

The entropy of the injected pairs is likely to increase over the initial entropy. The funnel 
is exposed to "acoustic" radiation from the accretion disk. A small fraction of the MHD 
waves generated by turbulence in the disk will propagate toward the funnel, be transmitted 
at the funnel wall, and dissipate within the funnel. Because the density of the funnel is so 
low, even a small fraction of the disk acoustic luminosity is capable of raising the pair plasma 
temperature to Be ^ 1 before it can escape at Lorentz factor Tjet ^ 1- Until the magnitude 



'^e± can be obtained by transforming from coordinate to fluid frame and subtracting the rest mass 
energy. 
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of turbulent heating can be estimated, the dynamics of the funnel pair plasma is extremely 
uncertain. 



5.4. Comparison to Goldreich-Julian density 

Pairs are created and then are accreted or escape on the light-crossing time T. This 
implies a density n± h±T in the funnel. Is n± sufficient to enforce the ideal MHD condition 
E = in the rest frame of the plasma {u^F^^ = 0): does the pair number density exceed 
the Goldreich-Julian density n^^j? 

A naive estimate of ncj uses the flat-space Goldreich-Julian charge number density 

^^■^ 47rec 3271^^/6 ^ ^ 

where the field rotation frequency in the funnel Q = (a*/8)c'^/ (GM) in the Blandford-Znajek 
model at a* < 1. For our standard Sgr A* model with B ~ 30G, a* ~ 0.94, M ~ 4.5x lO'^M©, 
so uqj ^ 10~'^cm~'^. 

A better estimate uses the Blandford-Znajek model for a monopole magnetosphere. We 
use Kerr-Schild coordinates (t, r, 9, cj)) and define the Goldreich-Julian density as the charge 
density measured in the frame of the normal observer, who to lowest (zeroth) order in a* 
has four-velocity = (—(1 + 2/r)~^ /^, 0, 0, 0). Using the c u rrent density J'^ derived from 



the BZ monopole solution as given in iMcKinney fc Gammid (120041 ) we find, to lowest order 
in a*, 

a,B'-c^ {l + 2/xy/^cos9 

" -"^^'^ = A^e ^''^ 

where x = r/C and B'' is the radial component of the field at z = 1. This is close to the 
naive estimate. @ Notice that h'qj ~ 1/r^, so n±/nGj is smallest closest to the horizon if 
n± ~ 1/r^ as, for example, in a wind. 

A still better estimate for ugj would use the simulation-derived currents. We have 
checked these and the charge densities are consistent with the estimates just given. 

Using estimates for B^ from § HI we can derive a scaling with rh and m: 

ncj ~ 3 X 10"^ a* m^/^ m" cm'^ (45) 



Since n'^j ~ cos9 the charge density changes sign from one hemisphere to the other. In a spht 
monopole model the sign would be the same. 
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near x = 2. Then using n± ~ n±T ~ 5 x lO^^m^ (using Equation [26] and /i^ = 1): 

^GJ n in-42 • -11/2 -3/2 //icN 

~ 6 X 10 a^, m ' mo ' . (46) 

n± 

This is the ratio at the axis. Although the Goldreich-Juhan density varies only weakly across 
the funnel, n±(/i^ = \^f)ln±{\x^ = 1) ^ 20 at x ~ 2. At very low m, where uqj ~ n±, the 
center of the funnel is populated by some other process — perhaps a pair cascade — and the 
edges by pair production in 77 collisions. 



6. Discussion 

6.1. Selfconsistency of the models 

Our models are self-consistent when they are radiatively inefficient and n± is greater 
than the Goldreich- Julian density. Figure [H] shows the self-consistent m, M as a shaded 
region. The region is bounded at low accretion rates by the solid line, where n± = ugj (for 
a^: = 0.94). The region is bounded at high accretion rates by rh = rhcru (vertical dashed line), 
where the model becomes radiatively efficient (here defined as Lboi/Mc^ = 0.1). Numerically, 
TTicrit ~ 10~^ at Ti/Te = 1. The models are fully self-consistent in the resulting wedge in 
parameter space. They are never applicable to stellar mass black holes, which cannot produce 
enough pairs to exceed the Goldreich- Julian density even at rhcrit- 



6.2. Sgr A* 

The bright radio source associated with the M = 4.5 x lO^Mo black hole in the Galactic 
center, Sgr A*, is a weak X-ray source ('quiescent' emission Ly < lO^^ergss"^) and is 



strongly sub-Eddington {L/LEdd ~ 10 ^). iMoscibrodzka et al.l ( 120091 ) have presented models 



of Sgr A* that suggest the most probable spin of the black hole a* = 0.94, temperature ratio 
Ti/Te = 3, and m = 2 x 10~^. The n± rate for these parameters is lower than the scaling laws 
of §5 because Ti/Tc > 1. For the purposes of this subsection only we consider 3D GRMHD 
models that have n± very similar to the 2D models. All models assume the accretion flow 
lies in the equatorial plane of the black hole. 

During the quiescent state the X-ray luminosity is Zx < 1- Near the horizon in the 
funnel, in model I, n± ^ 10~^s~^. The light crossing time is T ~ 20s, so a typical pair 
density near the horizon in the funnel is n± ^ 10~'^cm~^. This is flve orders of magnitude 
below Ugj ~ 10~^cm~^. In quiescence the funnel must therefore (for the assumed spin) be 
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populated by a process other than 77 pair production process considered here, for example 
a pair cascade. 

Sgr A* exhibits intraday variability at all observed wavelengths (radio, sub-mm, NIR, 
and X-rays). In particular in 2-10 keV luminosity r nay increase up to 16 times (the brightest 



flare detected has luminosity of Lx = 5.4 x 10*^^ iPorquet et al.ll2008l ) from the 'quiescent' 
level and last a few ks. During a bright flare the X-ray slope can change, and the implied 
pair density may reach or exceed ncj. Even during a flare the funnel kinetic luminosity is far 
below the Blandford-Znajek luminosity (see Table 1) for any reasonable Fjet. We conclude 
that close to the black hole any jet in Sgr A* is electromagnetically dominated. 



6.3. M87 



^The core of M87 hosts a sub-Eddingt on black hole with M = 3 x IO^Mq (IMarconi et al. 

19971 . but see iGebhardt fc ThomasI |2009| ) at distance D 16Mp c. M87 has a prominent 



radio jet resolved from lOOGM/c^ — Ikpc. [Reynolds et al.l (jl996bl ) have argued that models 
in which the jet is made of a pair plasma are favored over those in which the jet is composed 
of an ion-electron plasma. It is difficult to apply our model to M87 because the SED of M87 
from r < 100 GM/c^ includes contributions from both the accretion flow and the jet. 



We assume a*=0.94, set the inclination i = 30deg (IHeinz fc BegelmanI 119971 ) . as- 
sume that the accretion disk lies in the equatorial plane of the black hole, that the elec- 
tron distribution function is thermal, an d that T-jTp, = 3. The SED is normalized via 
f^{u = 230 GHz) = 1770 mJy at 16 Mpc jTan et al.lboosh . We flnd that the resulting zero 
cooling, 2D model with m = 1.5 x 10"^ (model K) is radiatively efficient (see Table 1) and 
therefore not self- consistent. 

To flnd a more self-consistent model we have run a GRMHD simulation with synchrotron 
and bremsstrahlung cooling, but not Compton cooling, included (model L, with m = 10~^). 
The cooling rates and cooling algorithm are presented in the Appendix. The efficiency is re- 
duced to ~ 30%. Figure [To] shows the SED for model L; the dotted line is the bremsstrahlung 
contribution, which is negligible. The model has Lx = lO^^ergss"^, L± 
and Lbz ^ lO^^ergss"^. 



7 X lO^^ergss ^, 



If we identify Lbz as the jet luminosity then the model is inconsistent with existing 

(I20091)'), which range 



estimates (see the us eful compilation of estimates in Table 3 of iLi et al 



from 3 X lO^^ergss ^ flYoung et al.ll2002l ) to > lO^^ergss ^ estimated by lBicknell fc Begelman 
(119961 ). The discrepancy between Lbz in model L and observations is by 1-3 orders of 



magnitude, but the lowest 'observed' value of Lbz could be possibly reached in a model 



- 19 - 



which combines T\/Tc > 1 and radiative coohng. 



The jet is optically thin to pair annhilation: t± ~ n±C(jT ~ 10 . It is also optically 
thick to pair production for TeV photons, T.y^{r) ~ aTTijjiC ~ 10^, (ri/ij ~ lO^^cm"^ is the 



infrared photon density calculated from the Monte Carlo simulations) 



The shape of the spatial distribution of pair production in model L is similar to that in 
models without cooling (although the scaling of the distribution changes). The implied pair 
density n± = h±T ~ 10 cm""^ (T ~ 10^ s), which is 10^ times larger than hgj 
in almost the entire computational domain. 



10- 



cm 



Because of the shortcomings of the model, however, it is useful to use a more nearly 
model- independent estimate of t he total pair pro d uction rate based on Equation fl^ . For 
Lx ~ 3 X 10"^^ (7 X 10^° from iDi Matteo et al.l (120031 ) corrected upward to an isotropic 
X-ray luminosity because our models beam X-rays into the equatorial plane) and ax = 0, 
A''-!- ~ lO^^s"^ This implies L/^ = fjetN±meC^Tjet = SxlO^^Tjetfjet- The implied pair density 
exceeds uqj for model L by ~ 10*. Since ugj oc -B oc rh^/^ and n± oc Lg^g the implied pair 
density will fall below the Goldreich- Julian density only for (m/10"^)^/^(L5i2/10^^-^)~^ < 10*. 
Even if m ~ 10~^ this would r equire L512 ^ 10*^*, wh ich seems implausibly low given the 



10"^" ergs s ^ TeV luminosity (lAharonian et al.ll2006l ). Therefore the main conclusion of 



this section does not change even if a more selfconsistent model is found. 

There are significant limitations on the model. We have considered only one value of 
a*; estimates and preliminary models not described here show that the pair production rate 
is a steeply increasing function of a^. Further preliminary models and a comparison of the 
Ti/Te = 3 model for Sgr A* with the scaling relation for Ti/Te = 1 models also show that the 
pair production rate declines sharply as Tj/Tp increases. But th e allowed values of Ti/Te are 
strongly constrained by submm VLBI (IFish &: Doelemarul2010l ). because as T^/T^ increases 
so does the size of the synchrotron photosphere. 



After submission of this article iLevinson fc Riegerl (120101 ) released a paper focused on 
modeling TeV emission and pair production in M87 (and Sgr A*). These authors use an 
ADAF model, assume that Tg saturates at few x lO^K (Og ~ 1), and set m ^ 10~^. The 
model is semi-analytic and does not include general relativistic effects. Bremsstrahlung is the 
dominant source of photons near the pair-production t hreshold, and the re s ulting radiation 
field is inadequate to raise the pair density above ugj. ILevinson fc Riegerl (120 lOl ) therefore 
invoke a gap/pair cascade model to produce pairs. 



We have investigated the ILevinson fc Riegerl (|2010l ) model by calculating images and an 
SED for a GRMHD/radiative transfer model with 0^ = 1 everywhere, m = 10~^, and a* = 
0.94. The model includes synchrotron, Compton and bremsstrahlung. We find f 230GHz = ^Jy 
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(at i = 30deg), and Lbz = lO^^ergss"^, consistent with observations. Free- free cooling 
dominates over synchrotron coohng only at r > 20 GM/c^. Levinson & Rieger neglect 
Compton cooling, but we find that Compton y = At ^ 12 and that with Compton cooling 
included the model efficiency is ~ 200%. The parameter space is large and the spectrum is 
parameter-sensitive, so there may be nearby models (with different a*, 0e, m, i) that are 
radiatively inefficient. The main point, however, is that self-consistent models can contain 
surprises that might not be anticipated in quasi-analytic estimates. Comptonization, in 
particular, occurs close to the innermost stable circular orbit, is therefore sensitive to the 
spin, and requires proper treatment of gravitational lensing. We concur with Levinson & 
Rieger's conclusion that the pair production rate due to 77 collisions is small. 

The model is also constrained by VLBI measurements. An optically thick spherical 
source of radius r and distance D in the Rayleigh- Jeans regime has ff ux ^ 27rQ pjUpC^ jr / D)" ^ / A^. 



Small r inferred from VLBI therefore requires high 6e. At 230 GHz (IFish &: Doelemanll2010l ) 
report structure on scales of a "few Schwarzschild" radii, while we find the Levinson & Rieger 
model has a photosphere at ~ SOGM/c^. In comparison, our model L has a photosphere 
at ~ IGM/c^. This argues against the Levinson & Rieger model if the reported structure 
arises from the accretion fiow rather than the jet. 



7. Summary 

We have studied electron-positron pair production in black hole magnetospheres by 77 
collisions. Our pair production rate simulations are based on a GRMHD time dependent 
model of a magnetized disk around a spinning black hole. The disk is a source of high 
energy radiation formed in multiple Compton scatterings of synchrotron photons. The pair 
production rates are calculated nearly ab-initio within 40GM/ of the event horizon, using 
Monte Carlo methods. 

The main results of this work are the fitting formulae for the rate and spatial distribution 
of pair production in terms of mg and rh (Equation [2H]) and in terms of mg, Lx, and a 
(Equation [30]). These indicate that 77 pair production is concentrated close to the event 
horizon, and is sensitive to model parameters such as rh. The pair production rate is also 
sensitive to black hole spin a* and the electron-ion temperature ratio Tj/Te, but exploring 
the dependence on these parameters is beyond the scope of this paper. 

We also find that the pair plasma is created with a power-law-like energy distribution. 
Most of the pairs are created in the equatorial plane of the thick disk because MeV photons 
created by Compton scattering are beamed into the equatorial plane. The pair plasma has 
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negligible_effect o n the accretion flow dynamic al evolution, consistent with previous results 
by lEsinlll999l and lKusunose fc Mineshigelll996l . assuming that it escapes on the viscous time 
scale. 

Only a few percent of all pairs are created in the magnetized funnel (black hole magne- 
tosphere), and most of pairs in the funnel are created near its wall. Pair jets will have spectra 
with a turnover frequency at around Vt = 10~^n-|-£ Hz (for example, for M87 £ = 4 x 10^^, 
and n± = 10, turnover frequency Ut = 10^^ Hz). 

We also flnd that the general relativistic RIAF models are self consistent up to rhr rit ~ 
10~^, which is consistent with the rhcrit = 5 x 10~® reported by [Fragile fc Meierl |2009| . For 
higher m one must couple the radiative cooling and forces into the dynamical model. 

Models with m < rhcrit have force-free, Thomson thin jets with the Blandford-Znajek 
luminosity much larger than pair kinetic luminosity. In models with very small m, the pair 
plasma density in the funnel is below the Goldreich- Julian density ncj, suggesting that 
another process, such as a pair cascades, will operate and populate the funnel. 

We have applied versions of our model to Sgr A* and to M87. These models suggest 
that n± > no J in M87 and n± < ncj in Sgr A*, with the important caveat that there 
are parameters (a* and Ti/Tc) that we have not varied, and effects (Compton cooling, and 
nonthermal electrons) that we have not included. 
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A. Synchrotron cooling rates including radiation self-absorption 

The synchrotron cooling rate for a single electron is 

,T ^ 2e^i?^(7^-l)sin^e 
3mgC^ 



This preprint was prepared with the A AS lATJrjX macros v5.0. 
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wher e ^ is the pitch angle between electron velocity and magnetic field (e.g. iRybicki fc Lightman 



19861 ). To obtain the total cooling rate from the thermal population of electrons we integrate 



Equation flAip against the relativistic Maxwellian distribution: 



--^^(^'-l)''%xp(--^) (A2) 



d-fd cos^ 2QJ<2{l/Qe) O 
The resulting integral over cos ^ and 7 is: 



e 



_ AB^eWej<,{l/e,) 

- Sc^mim/Q^) ^^^^ 



For X < 1, Kn{x) r(ra)/2(2/a;)'", so for large Ge 



This agree with expression 14 in IWardzihski fc Zdziarskil (j2000[ ). For a: ^ 1 (0e ^ 1), 

= Q 3 2 (^^) 
The ratio of these two expressions is 49e; a reasonable approximation is 

^ 4la/ti ^ + (26.)^-)^/- (A6) 

where m = 4/3 gives at most 4% error. 

To account for synchrotron self absorption, is multiplied by a factor: 

2 roc rTT 2 /"oo /•tt 

/ = - / du sm9d6j^exp{-T{u,9)) ^ / di^ sm9d9j^ (A7) 

^Syn Jo Jo ^Syn J v^rit Jo 

where is given by Equation (|3]), Asyn is the first integral without optical depth factor, 
Ucrit is the frequency where selfabsorption becomes important. The critical frequency is 
calculated numerically from: 

K^{9 = Tx/2)R = 1 (A8) 

where Ki, = j^/B^, By is the Planck function and R = O.IC. We find that / well approxi- 
mated by 

/4(exp(-^)+exp(-||i)) (A9) 

where Xcrit = ^cnt/i^s- This fit gives error for / less than 1% up to Xcnt = 10^ and 5% error 
at X,rit = 10^. 
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B. Free-free cooling 



The electron- ion bremsstrahlung cooling rate is fIStepney fc Guilbertlll983l ): 



^2 J f-(ln(2eeexp(-7i5) + 0.42) + 1.5) Qe > 1; 



4 



•20e \0.5/ 



1.7891-34) 



< Ge < 1. 



(Bl) 



where 7^ = 0.5772 is Euler constant a nd ar is the fin e structure constant. The electron- 
electron bremsstrahlung cooling rate is (jSvenssoru 119821 ) 



A 



- nV.ca.m I ^0e(ln(2e. exp(-7^)) + |) 
ee n^aTcajm^c <; ^(44 _ 3^2)01.5^1 + ^ ^q^ + @2 _ 1256^) 



ee>l; 

< 6e < 1. 



(B2) 

The cooling rates are i n units of erg ss'^ cm"^, and_are consistent within a factor of 2 with 
those provided by e.g. iMaxonI (119721 ) or iGouldl (jl980l ). Selfabsorption for free-free emission 
is negligible. For 6e > 1 the ratio of synchrotron to bremsstrahlung cooling rate is approxi- 
mately 6g//3a/. Synchrotron cooling dominates over the free-free emission in all of models 
considered here. 



C. Radiative cooling in MHD code 

Radiative cooling is governed by 

du du 1 A ^^^^ 
dt dr M* M* 

where u is the internal energy per unit proper volume, r is the proper time, and is the 
time component of the fiuid four-velocity. 

Numerically u is evolved in an operator-split fashion. After each fiuid timestep At, u is 
evolved using the second order scheme = u„exp(— At/rcooi,n+i/2) and TcooI = u/A. 

The cooling rates are calculated in cgs units, and then Acode = AcgsCT^ / M.. 



D. Bremsstrahlung emissivity in the radiative transfer calculations 



The emissivity for e-i interactions is (jStepney fc Guilbertlll983l ): 



dE 



3V 



dtdVdi/d^l 47r 



— rijch 



l+U} 



uj^l3nei'y)d'y 



(Dl) 
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where u = hu/meC^, ne{'~f) is relativistic Maxwellian electro n energy distribution an d the 
cross-section for this reactions is in the ultra-relativistic hmit fjjauch fc Rohrhchlll976l ). The 
l/Air factor gives emissivity per unit sohd angle. The integral is computed numerically using 
Gauss quadratures. The integration of Equation flDip over photon energies and solid angle 
gives the total e-i cooling rate, Agj. 



The emissivity for e-e emission is also from IStepney &: GuilbertI (119831 ) , 

1 



■ee 



—nlaTchajQe exp(— 6e) 



(D2) 



where x = {hiy/meC^)/Qe and G{x, 6e) is given in lStepney &: GuilbertI (119831 ). This formula 
is accurate to 5% over 0.1 < Og < 2. 



For Gg < 0.1 we use a quadrupole approximation ( iMaxonI 119721 ): 

1 2 



JT = 1 n^aTchafB[x] 



^exp(-f)J^o(f) 



(D3) 



where B{x) = 0.85 + 1.35^^ + 0.38a; and Kq is the modified Bessel function of the second 
kind. 



For 6e > 2 we use the ultra-relativistic approximation ( 1 AlexanianI 1 1 9 68l lMaxorull972[ ): 

(D4) 
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jr = -^t;^2(^Tcna f exp{-x){^ + 2x + ^ + 2(| + |x + x^) 
^'^^'^^ - 0.577] -exp(x)Ez(-a;)(| -43x ' 



x[lg( 



x')} 



Formulas ID2t ID3| ID4l connect smoothly at 0e=O.l and 2. Bremsstrahlung for e-e interactions 
dominates over e-i ones for Ge > 1 . Integration of j^^ over frequencies and solid angle gives 
the total cooling rate, Age. 



For details of the radiative transfer scheme see iDolence et al.l (120091 ): we sample the 
bremsstrahlung radiation field in the same way as for synchrotron radiation, except that 
bremsstrahlung is emitted isotropically in the fluid frame. For the range of parameters con- 
sidered in this work, energy loss by free-free emission is small in comparison to synchrotron 
and Compton losses. 
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Table 1. List of GRMHD models. 



ID 


a. 






< m >t 


LboI I Edd 


radiative 
efficiency 




note 


A 


0.94 


4.5 X 10" 


-2 


2 X 10-3 


10-11 


7 X 10-4 


10-17 


2D 


B 


0.94 


4.5 X 10" 


-2 


6 X 10-3 


10-10 


2 X 10--* 


10-15 


2D 


C 


0.94 


4.5 X 10" 


-2 


1 X 10-* 


4 X 10-10 


4 X 10--* 


10-13 


2D 


D 


0.94 


4.5 X 10" 


-2 


5 X 10-* 


2 X 10-* 


0.02 


10-10 


2D 


E 


0.94 


4.5 X 10' 


-2 


1 X 10-7 


6 X 10-* 


0.04 


10-* 


2D 


F 


0.94 


4.5 X 10" 


-3 


1 X 10-* 


5 X 10-10 


5 X 10--^ 


10-14 


2D 


G 


0.94 


4.5 X 10" 


-1 


1 X 10-* 


4 X 10-10 


4 X 10--' 


10-12 


2D 


H 


0.94 


4.5 




1 X 10-* 


3 X 10-* 


3 X 10--' 


10-11 


2D 












Sgr A* 








I 


0.94 


4.5 X 10" 


-2 


2.7 X 10-* 


5 X 10-10 


2 X 10-3 


10-11 


3D-quiescent,Ti/Tc = 3 


J 


0.94 


4.5 X 10" 


-2 


5.3 X 10-* 


1 X 10-3 


3 X 10-3 


10-3 


3D- weak flare, Ti/Tc = 3 












M87 








K 


0.94 


30 




1.5 X 10-'' 


3 X 10-* 


16.5 


0.1 


2D - w/o cooling 


L 


0.94 


30 




1 X lO-*' 


3 X lO-'^ 


0.3 


4 X 10-5 


2D - w/ cooling 



Note. — From left to right columns are: model ID, dimensionless spin of the black 
hole, the black hole mass in units of Mq, the rest mass accretion rate through 
the black hole horizon in units of Eddington mass accretion rate (MEdd = 2.22m8 
M0yr~^) averaged over later times of the simulation (At = 1500 — 2000T), the 
Eddington ratio LBoi/LEdd, the model radiative efficiency r] = Lboi/Mc^ {LboI 
is the RIAF luminosity integrated over emitting angles and frequencies), ratio of 
Kinetic to electromagnetic luminosity, and comments on models. Models I & J 
correspond to Sgr A* while K&L model M87. Run L accounts for cooling terms 
in the dynamical solution so the pair production rate is reduced. 
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Fig. 1. — Structure of RIAF. Panels from left to right: density distribution, plasma (3 
parameter, and dimensionless electron temperature, respectively. 
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Fig. 2. — Test problem: the pair production rate in the plane of two, isotropic point sources 
of high energy radiation {k\ = = ArrieC^). Pair production rate is given here in units of 
N^T/L^, where L is a length unit, A'^^ is a number of photons produced by each source per 
unit time T. Pair production rate is zero in two side regions because the energy of photons 
in the center-of-momentum frame is below the threshold energy there. The pair production 
rate is symmetric with respect to the axis connecting two sources. Black contour - see Fig|3l 




Fig. 3. — Test problem: upper panel: Analytical (green line) and numerical (red points) pair 
production rates, along the black contour in Figure El Lower panel: The fractional difference 
between analytical and numerical solutions. 
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Fig. 4. — Test problem: the grid averaged fractional difference between the numerical and 
analytical pair production rates as a function of number of photons packets produced by each 
source Ng. The dashed line is proportional to N^^^'^. For A^^ = 10'' the average difference 
per zone is about 7%, for A''^ = 10^ it is on average less than 3%. 
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Fig. 5. — Spatial distribution of h± in RIAF model C (points) and the contours of corre- 
sponding fitting function given by Equation |26] (lines). The fractional difference between 
model and data in this case is < 40%. Black contours mark the black hole horizon and the 
funnel wall. 
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Fig. 6. — Pair production rate dependence on the model parameters. Comparison of the 
total pair production rate N± to the fitting formula for models with various mass accretion 
rates m (A-E, blue filled symbols), and black hole masses m (F-H, red open symbols). The 
(analytical) is given by Equation l28l 
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Fig. 7. — Pair production rate dependence on observable parameters. Comparison of the 
total pair production rate N± to the fitting formula for models with different X-ray luminosi- 
ties {iyL^)2-iokeVy X-ray spectral index a and masses. Crosses, open squares, filled circles 
correspond to different snapshots in models A, B, and C, respectively. Open circles mark 
time averaged data from models with different masses (F, G, and H). The iV-t (analytical) is 
given by Equation [3ll 
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log 7 



Fig. 8. — The energy distribution of e"^ pairs produced in the magnetized funnel where 
7 = 7e±[FF] is measured in the plasma frame at different radii (single time slice of model C). 
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Fig. 9. — Fitting formulas are fully selfconsistent with the model assumptions for m and 
rh within the shaded region. Two solid lines mark regions where the ncj equals to the pair 
density at the funnel axis and funnel wall. Dotted lines shows scaling law for the number of 
scatterings (n^c)- Sgr A* and M87 are marked as open circles. 
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Fig. 10. — Model L: time averaged spectral energy distribution. Two lines show model 
with rh = 10~® and Ti/Tc = 1. m is chosen to normalize to 1.7 Jy at 230 GHz. Thick-solid 
line corresponds to run in which a radiative cooling is taken into account as described in 
the Appendix. Dotted lines is a free-free process spectrum. Model, that does not accounts 
for any cooling in the MHD simulation, is marked as thin-solid line (also this model i s 



not shown in th e table). 
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